﻿//[]-----------------------------------------------------------------------------------------------------------------[]
// | 开源程序: 使用SCE-UA算法对新安江水文模型的参数进行优化 （C#语言）
// |-----------------------------------------------------------------------------------------------------------------
// | 参考文献：
// | [0]河海大学水文水资源学院 卢家波 2021 https://github.com/lujiabo98/XAJ-SCEUA
// | [1]段青云,SCEUA的原始Fortran代码,1992, https://shxy.hhu.edu.cn/2019/0904/c12296a195177/page.htm
// | [2]L. Shawn Matott改编的C++代码,2009, https://github.com/MESH-Model/MESH_Project_Baker_Creek/blob/7e0a7e588213916deb2b6c11589df0d132d9b310/Model/Ostrich/SCEUA.h
// | [3]Van Hoey S改编的Python代码,2011
// | [4]Mostapha Kalami Heris, Shuffled Complex Evolution in MATLAB (URL: https://yarpiz.com/80/ypea110-shuffled-complex-evolution), Yarpiz, 2015.
// | [5]Duan, Q.Y., Gupta, V.K. & Sorooshian, S. Shuffled complex evolution approach for effective and efficient global minimization. J Optim Theory Appl 76, 501–521 (1993). https://doi.org/10.1007/BF00939380.
// | [6]Duan, Q., Sorooshian, S., and Gupta, V. (1992), Effective and efficient global optimization for conceptual rainfall-runoff models, Water Resour. Res., 28( 4), 1015– 1031, https://doi.org/10.1029/91WR02985.
// | [7]Duan, Q., Sorooshian, S., & Gupta, V. K. (1994). Optimal use of the SCE-UA global optimization method for calibrating watershed models. Journal of Hydrology, 158(3-4), 265-284. https://doi.org/10.1016/0022-1694(94)90057-4.
// | [8]王书功. 水文模型参数估计方法及参数估计不确定性研究[M]. 河南:黄河水利出版社,2010.(https://book.douban.com/subject/5377630/)
// | [9]王书功. 水文模型参数估计方法及参数估计不确定性研究[D]. 北京:中国科学院研究生院,2006.(https://jz.docin.com/p-87849994.html)
//[]-----------------------------------------------------------------------------------------------------------------[]
/*
考虑最小化问题

待优化参数信息:
        NOPT = 待优化的参数数量，问题维数，目标函数含有的变量数，简记为n
        XNAME(.) = 变量名
		A(.) = 初始参数集
		BL(.) = 参数下限
		BU(.) = 参数上限
		
SCE-UA算法控制参数：
		NGS = 初始种群中的复形数量，简记为p
		NPG = 每个复形中的点的个数，简记为m		
		NPS = 子复形中的点数，简记为q
        ALPHA = CCE步骤4中 a-e 重复的次数，即复形洗牌前每个复形允许的进化步骤数
        BETA = CCE步骤3-5重复的次数	
        NPT = 初始种群中的总点数(NPT=NGS*NPG)，简记为s

收敛检查参数：
		MAXN = 优化终止前允许的最大试验次数
		KSTOP = 在终止优化之前，标准值必须改变给定百分比的洗牌循环数
		PCENTO = 给定数量的随机循环中，标准值必须改变的百分比         
		PEPS = 给定数量的随机循环中，函数值最低变化率

标志变量:
		INIFLG = 标记是否将初始点包括在种群中
			= false, 不包括
			= true, 包括
		IPRINT = 用于控制每个洗牌循环后的打印输出的标志
			= false, 打印关于种群最佳点的信息
			= true, 打印关于种群每个点的信息
        IDEFLT = 是否采用默认参数值标志
		    = false, 采用
			= true, 不采用

优化结果：
        BESTX(.,.) = 每次洗牌循环的最佳点
        BESTF(.) = 每次循环的最佳点函数值

收敛判据：
		ICALL = 模型调用次数
		TIMEOU = 最近KSTOP次数函数值变化率
		GNRNG = 参数范围的归一化几何平均值
*/
//[]-----------------------------------------------------------------------------------------------------------------[]

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Threading.Tasks;

namespace SceuaApp
{
    public class SCEUA
    {
        // 每个变量的参数名、初始值、下限、上限
        public int m_nopt;  // 待优化的参数数量，问题维数，目标函数含有的变量数，简记为n
        public List<string> m_xname = new List<string>(); // 变量名
        public List<double> m_a = new List<double>();   // 参数初始值
        public List<double> m_bl = new List<double>();  // 参数下限
        public List<double> m_bu = new List<double>();  // 参数上限

        // SCE控制参数
        public int m_ngs;   // 种群中的复形数量，简记为p
        public int m_npg;   // 每个复形中的点的个数，简记为m
        public int m_nps;   // 子复形即单纯形中的点数，简记为q
        public int m_alpha;   // CCE步骤4中 a-e 重复的次数，即复形洗牌前每个复形允许的进化步骤数
        public int m_beta;    // CCE步骤3-5重复的次数
        public int m_npt;   // 种群中的总点数(NPT=NGS*NPG)，简记为s

        // 收敛判据参数
        public int m_maxn;    // 优化终止前允许的最大试验次数
        public int m_kstop;   // 在终止优化之前，计算函数值变化率的洗牌循环数
        public double m_pcento;  // 给定数量的随机循环中，函数值最低变化率
        public double m_peps;    // 判断所有参数收敛的标准

        // 标志变量
        public bool m_iniflg;  // 标记是否将初始点包括在种群中
        public bool m_iprint;  // 用于控制每个洗牌循环后的打印输出的标志
        public int m_ideflt;  // 采用默认值标志

        // 存储结果
        public List<List<double>> m_bestx = new List<List<double>>(); // 每次洗牌循环的最佳点
        public List<double> m_bestf = new List<double>(); // 每次循环的最佳点函数值

        // 存储收敛判据
        public List<int> m_icall = new List<int>();  // 优化目标试验次数
        public List<double> m_timeou = new List<double>();  // 最近kstop次洗牌循环中，函数值变化率
        public List<double> m_gnrng = new List<double>();  // 参数范围的归一化几何平均值

        // 存储待率定模型实测值与模拟值
        public List<double> measuredValues = new List<double>();
        public List<double> simulatedValues = new List<double>();

        // 存储SCEUA程序命令行参数
        public int argc;  // 参数个数
        public List<string> argv = new List<string>();  // 参数名
        public string filePath; // 可执行程序所在目录


        /********************************************************
        Optimize()
        使用 SCEUA 最小化目标函数
        *********************************************************/
public void Optimize()
        {
            scemain(); // 执行SCE程序
        }

        /**********************************************************
        scemain()
        执行SCE程序，调用scein()、sceua()、sceout()函数
        ***********************************************************/
        public void scemain()
        {
            scein();   // 设置优化参数
            sceua();   // SCE-UA算法
            sceout();  // 输出优化结果
        }

        /************************************************************
        scein()
        该函数读取用于全局优化的SCE方法的输入变量

        输入参数列表：
            StrPath = 文件路径，默认为空，即可执行程序所在路径
        ************************************************************/
        public void scein()
        {
            // 初始化 I/O 变量
            string filePathWithFile = Path.Combine(filePath, "scein.txt");
            using (StreamReader fin = new StreamReader(filePathWithFile))
            {
                // 读取参数个数
                m_nopt = int.Parse(fin.ReadLine());

                // 创建参数名、初始值、上下限数组
                m_xname = new List<string>(new string[m_nopt]);
                m_a = new List<double>(new double[m_nopt]);
                m_bl = new List<double>(new double[m_nopt]);
                m_bu = new List<double>(new double[m_nopt]);

                // 依次读取每个变量的参数名、初始值、下限、上限
                for (int i = 0; i < m_nopt; i++)
                {
                    var line = fin.ReadLine().Split(new[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    m_xname[i] = line[0]; // 读取变量名
                    m_a[i] = double.Parse(line[1]); // 读取参数初始值
                    m_bl[i] = double.Parse(line[2]); // 读取参数下限
                    m_bu[i] = double.Parse(line[3]); // 读取参数上限
                }

                // 读取默认值标识
                m_ideflt = int.Parse(fin.ReadLine());
                //fin.ReadLine();

                // 如果 IDEFLT 等于 true，接着读入 SCE 控制参数、收敛判据参数、标志变量
                if (m_ideflt == 1)
                {
                    var line = fin.ReadLine().Split(new[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    m_ngs = int.Parse(line[0]);
                    m_npg = int.Parse(line[1]);
                    m_nps = int.Parse(line[2]);
                    m_alpha = int.Parse(line[3]);
                    m_beta = int.Parse(line[4]);

                    line = fin.ReadLine().Split(new[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    m_maxn = int.Parse(line[0]);
                    m_kstop = int.Parse(line[1]);
                    m_pcento = double.Parse(line[2]);
                    m_peps = double.Parse(line[3]);

                    line = fin.ReadLine().Split(new[] { ' ', '\t' }, StringSplitOptions.RemoveEmptyEntries);
                    m_iniflg = bool.Parse(line[0]);
                    m_iprint = bool.Parse(line[1]);
                }
                // 否则将 SCE 控制参数设置为推荐值即默认值
                else
                {
                    m_ngs = 5;
                    m_npg = 2 * m_nopt + 1;
                    m_nps = m_nopt + 1;
                    m_alpha = 1;
                    m_beta = 2 * m_nopt + 1;
                    m_maxn = 5000;
                    m_kstop = 10;
                    m_pcento = 0.1;
                    m_peps = 0.001;

                    m_iniflg = true;
                    m_iprint = false;
                }

                // 计算初始种群中的总点数
                m_npt = m_ngs * m_npg;
            }

            // 读取实测值
            measuredValues = ReadValues(Path.Combine(filePath, "observe.txt"));
        }

        /**********************************************************
        sceout()
        输出SCE-UA算法运行结果

        输入参数列表：
            StrPath = 文件路径，默认为空，即可执行程序所在路径

        局部变量列表：
            fout = 文件输出句柄
        ***********************************************************/
        public void sceout()
        {
            // ==================1. 控制台输出=================
            Console.WriteLine("===================参数初始值及上下限===================");
            Console.Write("参数名：");
            foreach (var name in m_xname)
            {
                Console.Write(name + "  ");
            }
            Console.WriteLine();

            Console.Write("参数初始值：");
            foreach (var value in m_a)
            {
                Console.Write(value + "  ");
            }
            Console.WriteLine();

            Console.Write("参数下限：");
            foreach (var lower in m_bl)
            {
                Console.Write(lower + "  ");
            }
            Console.WriteLine();

            Console.Write("参数上限：");
            foreach (var upper in m_bu)
            {
                Console.Write(upper + "  ");
            }
            Console.WriteLine();

            Console.WriteLine("初始参数函数值：" + functn(m_nopt, m_a));

            Console.WriteLine("===================SCE-UA搜索的结果===================");
            Console.WriteLine("最优点:");

            foreach (var name in m_xname)
            {
                Console.Write(name + "  ");
            }
            Console.WriteLine();

            //在C#中，可以使用 List<T>.Last() 方法来获取 List<T> 的最后一个元素。
            //需要注意的是，Last() 是 LINQ 方法，因此需要在文件顶部添加 using System.Linq;
            foreach (var value in m_bestx.Last())
            {
                Console.Write(value + "  ");
            }
            Console.WriteLine();

            Console.WriteLine("最优函数值:" + m_bestf.Last());
            Console.WriteLine(m_bestf.Count + "次洗牌演化的最优点函数值：");

            foreach (var value in m_bestf)
            {
                Console.Write(value + "  ");
            }
            Console.WriteLine();

            // ==================2. 文件输出=================
            string filePathWithFile = Path.Combine(filePath, "sceout.txt");
            using (StreamWriter fout = new StreamWriter(filePathWithFile))
            {
                fout.WriteLine("===================参数初始值及上下限===================");
                fout.Write("参数名：");
                foreach (var name in m_xname)
                {
                    fout.Write(name + "  ");
                }
                fout.WriteLine();

                fout.Write("参数初始值：");
                foreach (var value in m_a)
                {
                    fout.Write(value + "  ");
                }
                fout.WriteLine();

                fout.Write("参数下限：");
                foreach (var lower in m_bl)
                {
                    fout.Write(lower + "  ");
                }
                fout.WriteLine();

                fout.Write("参数上限：");
                foreach (var upper in m_bu)
                {
                    fout.Write(upper + "  ");
                }
                fout.WriteLine();

                fout.WriteLine("初始参数函数值：" + functn(m_nopt, m_a));

                fout.WriteLine("===================SCE-UA搜索的结果===================");
                fout.WriteLine("最优点:");

                foreach (var name in m_xname)
                {
                    fout.Write(name + "  ");
                }
                fout.WriteLine();

                foreach (var value in m_bestx.Last())
                {
                    fout.Write(value + "  ");
                }
                fout.WriteLine();

                fout.WriteLine("最优函数值:" + m_bestf.Last());
                fout.WriteLine(m_bestf.Count + "次洗牌演化的最优点函数值：");

                foreach (var value in m_bestf)
                {
                    fout.Write(value + "  ");
                }
                fout.WriteLine();
            }
        }



        /*************************************************************
         sceua()
         该算法对目标函数最小化问题进行全局优化搜索

         局部变量列表：
             X(.,.) = 种群中点的坐标
             XF(.) = X(.,.)的函数值
             CX(.,.) = 复形中点的坐标
             CF(.) = CX(.,.)的函数值
             XNSTD(.) = 种群中参数的标准差
             BOUND(.) = 优化变量的约束
             ICALL = 模型调用次数
             TIMEOU = 最近kstop次洗牌循环中函数值变化率
             GNRNG = 参数范围的归一化几何平均值
             NLOOP = 主循环次数
        *************************************************************/
        public void sceua()
        {
            Console.WriteLine("==================================================");
            Console.WriteLine("                  进入SCE-UA全局搜索                ");
            Console.WriteLine("==================================================");

            //===========1.[初始化]定义局部变量及数组，并分配内存===========
            List<List<double>> x = new List<List<double>>(); // 种群中点的坐标
            for (int i = 0; i < m_npt; i++)
            {
                x.Add(new List<double>(new double[m_nopt]));
            }

            List<double> xf = new List<double>(new double[m_npt]); // 种群中点的坐标X(.,.)的函数值

            List<List<double>> cx = new List<List<double>>(); // 复形中点的坐标
            for (int i = 0; i < m_npg; i++)
            {
                cx.Add(new List<double>(new double[m_nopt]));
            }

            List<double> cf = new List<double>(new double[m_npg]); // 复形中点的坐标CX(.,.)的函数值

            List<double> xnstd = new List<double>(new double[m_nopt]); // 种群中参数的标准差，初始化为0.1
            for (int i = 0; i < m_nopt; i++)
            {
                xnstd[i] = 0.1;
            }

            List<double> bound = new List<double>(new double[m_nopt]); // 优化变量的约束范围

            int icall = 0; // 模型调用次数
            double timeou = 10000.0; // 最近kstop次洗牌循环中函数值变化率，默认为一较大数值
            double gnrng = 10000.0; // 参数范围的归一化几何平均值，默认为一较大数值
            int nloop = 0; // 主循环次数

            for (int j = 0; j < m_nopt; j++)
            {
                bound[j] = m_bu[j] - m_bl[j]; // 计算待优化参数约束范围
            }

            //===========2.[生成样本]在参数空间内生成一个含有NPT个点的初始集合===========
            GenerateSample(x, xf, ref icall);

            //===========3.[样本点排序]样本点按函数值升序排列===========
            RankPoints(x, xf);

            Console.WriteLine("目标函数调用次数" + "  " + "函数值变化率" + "  " + "参数变化范围");

            object lockObject = new object(); // 用于同步的锁对象

            //===========7.[收敛判断]不满足收敛判据，返回步骤4===========
            while ((icall < m_maxn) && (timeou > m_pcento) && (gnrng > m_peps))
            {
                nloop += 1; // 主循环次数计数

                // 对每个复形进行独立演化计算，并放回到缓冲区D中  
                //for (int k = 0; k < m_ngs; k++)
                Parallel.For(0, m_ngs, k =>
                {
                    lock (lockObject)
                    {
                        //===========4.[划分复形群体]将样本点划分到p个复形中，每个复形含m个点===========
                        Partition2Complexes(k, x, xf, cx, cf);

                        //===========5.[复形演化]根据CCE对每个复形独立演化计算===========
                        cce(cx, cf, xnstd, ref icall);

                        //===========6.[复形洗牌]将所有复形内的点放回缓冲区D===========
                        ShuffleComplexes(k, x, xf, cx, cf);
                    }
                });

                //===========3.[样本点排序]将所有复形内的点按函数值升序排列===========
                RankPoints(x, xf);

                //===========7.[收敛判断]如果满足收敛判据，则算法终止，否则回到步骤4===========
                CheckConvergence(x, xnstd, bound, nloop, icall, ref timeou, ref gnrng);

                Console.WriteLine("      " + icall + "            " + timeou + "        " + gnrng);
            }
        }


        /************************************************************
    GenerateSample()
    2.[生成样本]在参数空间内生成一个含有NPT个点的初始集合

    输入参数列表：   
        X(.,.) = 种群中点的坐标
        XF(.) = X(.,.)的函数值
        ICALL = 模型调用次数

    局部变量列表：
        XX(.) = 种群中点X(.,.)中单点的坐标
    ************************************************************/
        public void GenerateSample(
            List<List<double>> x,
            List<double> xf,
            ref int icall
        )
        {
            List<double> xx = new List<double>(new double[m_nopt]); // 种群中点X(.,.)中单点的坐标

            // 生成在参数空间均匀分布的NPT个随机点
            for (int i = 0; i < m_npt; i++)
            {
                // 在参数空间内按照均匀分布随机生成点
                getpnt(xx);

                // 将生成的点赋值给种群
                x[i] = new List<double>(xx);
            }

            // 如果INIFLG等于true，设置X(0,.)为初始点A(.)，将初始点包含在初始集合中
            if (m_iniflg == true)
            {
                x[0] = new List<double>(m_a);
            }

            // 计算NPT个随机点对应的函数值
            for (int i = 0; i < m_npt; i++)
            {
                xf[i] = functn(m_nopt, x[i]); // 计算生成点的函数值
                icall += 1; // 模型调用次数增加1
            }
        }

        /************************************************************
        RankPoints()
        3.[样本点排序]样本点按函数值升序排列

        输入参数列表：
            X(.,.) = 种群中点的坐标
            XF(.) = X(.,.)的函数值
        ************************************************************/
        public void RankPoints(
            List<List<double>> x,
            List<double> xf
        )
        {
            sort(x, xf); // 根据目标函数值将样本点按从小到大排序

            // 记录最好的点
            m_bestx.Add(new List<double>(x[0])); // 最好点的坐标，即第一个点
            m_bestf.Add(xf[0]); // 最好点的函数值，即第一个点的函数值
        }

        /************************************************************
        Partition2Complexes()
        4.[划分复形群体]将样本点划分到ngs个复形中，每个复形含npg个点

        输入参数列表：
            K = 复形的序号，K = 1,2,···,ngs
            X(.,.) = 种群中点的坐标
            XF(.) = X(.,.)的函数值
            CX(.,.) = 复形中点的坐标
            CF(.) = 复形中点的坐标CX(.,.)的函数值

        局部变量列表：
            INDEX = 按照规则（k + ngs * j）从种群中提取点到复形中的下标
        ************************************************************/
        public void Partition2Complexes(
            int k,
            List<List<double>> x,
            List<double> xf,
            List<List<double>> cx,
            List<double> cf
        )
        {
            int index = 0; // 种群的下标

            // 遍历复形中的每个点，为其从种群中赋值
            for (int j = 0; j < m_npg; j++)
            {
                index = k + m_ngs * j; // 计算下标

                cx[j] = new List<double>(x[index]); // 划分样本点到复形
                cf[j] = xf[index]; // 划分对应的函数值
            }
        }


        /************************************************************
         cce()
         5.[复形演化]根据CCE对每个复形独立演化计算

         输入参数列表：
            CX(.,.) = 复形中点的坐标
            CF(.) = 复形中点的坐标CX(.,.)的函数值
            XNSTD(.) = 种群中参数的归一化标准差
            ICALL = 模型调用次数

         局部变量列表：
            s(.,.) = 当前单纯形中点的坐标
            sf(.) = 当前单纯形中点的坐标S(.,.)的函数值
            sb(.) = 单纯形的最佳点
            sw(.) = 单纯形的最坏点WO（WORST POINT，简称WO）
            ce(.) = 单纯形排除最坏点（WO）的形心
            snew(.) = 从单纯形生成的新点
            par(.) = 除最差点以外的参数取值集合
            lcs(.) = 索引定位S(.,.)在 CX(.,.)中的位置
            wts = 复形中每个点的权重,weights
            vals = 经过均匀随机扰动的复形中每个点的权重
            valsWithIndices = 权重及下标索引
            gen = 随机数引擎
            u = 产生[0,1)均匀分布随机数
            fw = 最差点的函数值
            fnew = 新点SNEW的函数值    
            ibound = 违反约束指标
                    = false,  没有违反
                    = true,  违反
        ************************************************************/
        public void cce(
            List<List<double>> cx,
            List<double> cf,
            List<double> xnstd,
            ref int icall
        )
        {
            //===========1.[初始化]定义局部变量===========

            // 定义数组并分配内存
            var s = new List<List<double>>(Enumerable.Repeat(new List<double>(new double[m_nopt]), m_nps).ToList()); // 当前单纯形中点的坐标
            var sf = new List<double>(new double[m_nps]); // 当前单纯形中点的坐标S(.,.)的函数值
            var sb = new List<double>(new double[m_nopt]); // 子复形（即单纯形）的最佳点，其函数值最小
            var sw = new List<double>(new double[m_nopt]); // 子复形（即单纯形）的最差点，其函数值最大
            var ce = new List<double>(new double[m_nopt]); // 子复形（即单纯形）排除最坏点（WO）的形心
            var snew = new List<double>(new double[m_nopt]); // 从子复形（即单纯形）生成的新点
            var par = new List<double>(new double[m_nps - 1]); // 除最差点以外的参数取值集合
            var lcs = new List<int>(new int[m_nps]); // 索引定位S(.,.)在 CX(.,.)中的位置
            var wts = new List<int>(new int[m_npg]); // 用于存储复形中每个点的权重
            var vals = new List<double>(); // 用于存储经过随机扰动的复形中每个点的权重
            var valsWithIndices = new List<(int Index, double Value)>(); // 存储下标索引及其权重

            Random gen = new Random(); // 随机数引擎
            Func<double> u = () => gen.NextDouble(); // 产生[0,1)均匀分布随机数

            double fw = 0.0; // 最差点的函数值
            double fnew = 0.0; // 新点SNEW的函数值
            bool ibound = true; // 违反约束指标，默认违反

            //===========2.[分配权重]根据三角概率分布，为复形内每个点分配概率===========
            for (int i = 0; i < m_npg; i++)
            {
                wts[i] = m_npg - i;
            }

            //===========6.[迭代]重复步骤 3-5 beta次===========
            for (int ibeta = 0; ibeta < m_beta; ibeta++)
            {
                //===========3.[选择父辈群体]从复形中加权不放回抽取q个点===========
                vals.Clear();
                valsWithIndices.Clear();

                foreach (var iter in wts)
                {
                    vals.Add(Math.Pow(u(), 1.0 / iter));
                }

                for (int i = 0; i < m_npg; i++)
                {
                    valsWithIndices.Add((i, vals[i]));
                }

                valsWithIndices = valsWithIndices.OrderByDescending(x => x.Value).ToList();

                for (int i = 0; i < m_nps; i++)
                {
                    lcs[i] = valsWithIndices[i].Index;
                }

                lcs.Sort();

                for (int i = 0; i < m_nps; i++)
                {
                    s[i] = new List<double>(cx[lcs[i]]);
                    sf[i] = cf[lcs[i]];
                }

                //===========4.[产生下一代群体]按照下山单纯形法的6步对q个点演化计算,重复alpha次===========
                for (int ialpha = 0; ialpha < m_alpha; ialpha++)
                {
                    sb = new List<double>(s.First());
                    sw = new List<double>(s.Last());
                    fw = sf.Last();

                    for (int j = 0; j < m_nopt; j++)
                    {
                        for (int i = 0; i < m_nps - 1; i++)
                        {
                            par[i] = s[i][j];
                        }

                        ce[j] = par.Average();
                    }

                    // 计算新点SNEW
                    for (int j = 0; j < m_nopt; j++)
                    {
                        snew[j] = 2 * ce[j] - sw[j];
                    }

                    chkcst(snew, ref ibound);

                    if (ibound == true)
                    {
                        getpnt(sb, xnstd, snew);
                    }

                    fnew = functn(m_nopt, snew);
                    icall += 1;

                    if (fnew < fw)
                    {
                        s[s.Count - 1] = new List<double>(snew);
                        sf[sf.Count - 1] = fnew;
                    }
                    else
                    {
                        for (int j = 0; j < m_nopt; j++)
                        {
                            snew[j] = (ce[j] + sw[j]) / 2.0;
                        }

                        fnew = functn(m_nopt, snew);
                        icall += 1;

                        if (fnew < fw)
                        {
                            s[s.Count - 1] = new List<double>(snew);
                            sf[sf.Count - 1] = fnew;
                        }
                        else
                        {
                            getpnt(sb, xnstd, snew);

                            fnew = functn(m_nopt, snew);
                            icall += 1;

                            s[s.Count - 1] = new List<double>(snew);
                            sf[sf.Count - 1] = fnew;
                        }
                    }

                    sort(s, sf);
                }

                //===========5.[以子代取代父代群体]将q个点按照位置放回复形，升序排列===========
                for (int i = 0; i < m_nps; i++)
                {
                    cx[lcs[i]] = new List<double>(s[i]);
                    cf[lcs[i]] = sf[i];
                }

                sort(cx, cf);
            }
        }


        /************************************************************
        ShuffleComplexes()
        6.[复形洗牌]将所有复形内的点放回缓冲区D

        输入参数列表：
            K = 复形的序号，K = 1,2,...,ngs
            X(.,.) = 种群中点的坐标
            XF(.) = X(.,.)的函数值
            CX(.,.) = 复形中点的坐标
            CF(.) = 复形中点的坐标CX(.,.)的函数值

        局部变量列表：
            INDEX = 按照规则（k + ngs * j）从复形中提取点到种群中的下标
        ************************************************************/
        public void ShuffleComplexes(
            int k,
            List<List<double>> x,
            List<double> xf,
            List<List<double>> cx,
            List<double> cf
        )
        {
            int index = 0; // 种群的下标

            // 遍历复形中的每个点，将其值赋给种群
            for (int j = 0; j < m_npg; j++)
            {
                index = k + m_ngs * j; // 计算下标

                x[index] = new List<double>(cx[j]); // 样本点从复形赋值到种群
                xf[index] = cf[j]; // 函数值从复形赋值到种群
            }
        }

        public void CheckConvergence(
            IReadOnlyList<IReadOnlyList<double>> x,  // 种群坐标
            IList<double> xnstd,                    // 归一化标准差（输出参数）
            IReadOnlyList<double> bound,             // 变量约束
            int nloop,                               // 主循环次数
            int icall,                               // 模型调用次数
            ref double timeou,                       // 函数值变化率（输出参数）
            ref double gnrng                         // 参数范围几何均值（输出参数）
        )
        {
            //==========1. 检查函数评估次数上限==========
            if (icall > m_maxn)
            {
                Console.WriteLine($"经过{nloop}次洗牌演化，" +
                    $"优化搜索已终止，超过最大试验次数{m_maxn}次限制");
            }

            //==========2. 检查函数值停滞情况==========
            double denomi = 0.0;

            if (nloop >= m_kstop)
            {
                // 计算最近kstop次目标函数绝对值的均值
                var recentBestF = m_bestf.Skip(m_bestf.Count - m_kstop).Take(m_kstop);
                denomi = Math.Abs(recentBestF.Sum()) / m_kstop;

                // 计算函数值变化率
                double currentBest = m_bestf[nloop - 1];
                double pastBest = m_bestf[nloop - m_kstop];
                timeou = Math.Abs(currentBest - pastBest) / denomi;

                if (timeou < m_pcento)
                {
                    Console.WriteLine($"最佳点在最近{m_kstop}次循环中" +
                        $"函数值变化率小于阈值{m_pcento}");
                    Console.WriteLine($"经过{nloop}次洗牌演化，" +
                        "基于目标函数标准已实现收敛！");
                }
            }

            //==========3. 检查种群聚集程度==========
            double[] xmax = new double[m_nopt];
            double[] xmin = new double[m_nopt];
            double[] xmean = new double[m_nopt];
            double[] par = new double[m_npt];
            const double delta = 1.0e-20;
            double gsum = 0.0;

            for (int k = 0; k < m_nopt; k++)
            {
                // 提取当前参数所有取值
                for (int i = 0; i < m_npt; i++)
                {
                    par[i] = x[i][k];
                }

                xmax[k] = par.Max();
                xmin[k] = par.Min();
                double parsum = par.Sum();
                xmean[k] = parsum / m_npt;

                // 计算方差和归一化标准差
                double parsum2 = par.Sum(val => val * val);
                xnstd[k] = Math.Sqrt(parsum2 / m_npt - Math.Pow(xmean[k], 2));
                xnstd[k] /= bound[k];

                // 累加对数范围
                gsum += Math.Log(delta + (xmax[k] - xmin[k]) / bound[k]);
            }

            // 计算几何平均值
            gnrng = Math.Exp(gsum / m_nopt);

            if (gnrng < m_peps)
            {
                Console.WriteLine($"经过{nloop}次洗牌演化，" +
                    "种群已收敛到指定小参数空间");
            }

            //==========4. 记录收敛指标==========
            m_icall.Add(icall);
            m_timeou.Add(timeou);
            m_gnrng.Add(gnrng);
        }

        /**********************************************************
         getpnt()
         在可行区域内根据均匀分布生成一个新点

         输入参数列表：
             SNEW(.) = 产生的新点
        ***********************************************************/
        public void getpnt(List<double> snew)
        {
            // 随机数引擎采用设备熵值保证随机性
            var gen = new Random();

            // 违反约束指标，默认违反
            bool ibound = true;

            while (ibound == true)
            {
                // 根据要求产生新点
                for (int j = 0; j < m_nopt; j++)
                {
                    // 产生均匀分布随机数
                    double u = gen.NextDouble() * (m_bu[j] - m_bl[j]) + m_bl[j];
                    snew[j] = u;
                }

                // 显性和隐性检查
                chkcst(snew, ref ibound);
            }
        }

        /**********************************************************
         getpnt()
         在可行区域内根据正态分布生成一个新点

         输入参数列表：
             SNEW(.) = 产生的新点
             std(.) = 概率分布的标准差
             xi(.) = 焦点
        ***********************************************************/
        public void getpnt(List<double> xi, List<double> std, List<double> snew)
        {
            // 随机数引擎采用设备熵值保证随机性
            var gen = new Random();

            // 显性和隐性检查
            bool ibound = true; // 违反约束指标，默认违反

            while (ibound == true)
            {
                // 根据要求产生新点
                for (int j = 0; j < m_nopt; j++)
                {
                    // 产生正态分布随机数
                    double n = NextNormal(gen, xi[j], std[j]);
                    snew[j] = n;
                }

                // 显性和隐性检查
                chkcst(snew, ref ibound);
            }
        }

        private double NextNormal(Random rng, double mean, double stdDev)
        {
            double u1 = rng.NextDouble();
            double u2 = rng.NextDouble();
            double randStdNormal = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
            return mean + stdDev * randStdNormal;
        }

        /************************************************************
         sort()
         按函数值的升序对一组点数组和对应的函数值数组进行排序

         输入参数列表：
             X(.,.) = 种群中所有点的坐标
             XF(.) = X(.,.)的函数值
        *************************************************************/
        public void sort(List<List<double>> x, List<double> xf)
        {
            // 将点的坐标与函数值对应上，以便于扩展排序
            var xxf = new List<(List<double>, double)>();

            // 构造向量xxf，存储坐标和函数值对
            for (int i = 0; i < xf.Count; i++)
            {
                xxf.Add((x[i], xf[i]));
            }

            // 匿名函数定义比较规则，用于pair数据结构按照函数值升序排列
            xxf.Sort((a, b) => a.Item2.CompareTo(b.Item2));

            // 将排序后的xxf赋值给点数组x和函数值数组xf
            for (int i = 0; i < xxf.Count; i++)
            {
                x[i] = xxf[i].Item1;  // 取出点坐标
                xf[i] = xxf[i].Item2; // 取出函数值
            }
        }

        /************************************************************
         chkcst()
         检查试验点是否满足所有约束

         输入参数列表：
             XX(.) = X 中单点的坐标
             IBOUND = 违反约束指标
                     = false,  没有违反
                     = true,  违反
        *************************************************************/
        public void chkcst(List<double> xx, ref bool ibound)
        {
            ibound = true; // 默认违反约束

            // 检查是否违反了显式边界约束
            for (int i = 0; i < m_nopt; i++)
            {
                if ((xx[i] < m_bl[i]) || (xx[i] > m_bu[i]))
                {
                    // 至少违反了一个约束
                    return; // 终止函数运行
                }
            }

            // 检查是否违反了隐式约束
            // 此处用于填写隐式约束
            // 此函数没有隐式约束

            // 通过上述检查，表明没有违反约束
            ibound = false;
        }

        /************************************************************
         functn()
         计算目标函数值

         输入参数列表：
             NOPT = 待优化的参数数量，问题维数，目标函数含有的变量数，简记为n
             X(.) = 单点的坐标，即各参数构成的数组
        *************************************************************/
        public double functn(int nopt, List<double> x)
        {
            // 1.[前处理]在functn中先把自动生成的参数写入到新安江模型的输入文件中
            PreProcessing(x);

            // 2.[运行模型]再调用新安江模型计算出结果
            RunModel();

            // 3.[后处理]最后计算出NSE，将1-NSE作为functn函数返回值
            double returnVal = PostProcessing();

            return returnVal;
        }

        public void PreProcessing(List<double> x)
        {
            // 初始化 I/O 变量
            string filePath = "parameter.tpl";
            string outputPath = "parameter.txt";

            using (var fin = new System.IO.StreamReader(filePath))
            using (var fout = new System.IO.StreamWriter(outputPath))
            {
                string parameter = "";
                while ((parameter = fin.ReadLine()) != null)
                {
                    // 查找是否在参数名数组中
                    int index = m_xname.IndexOf(parameter);
                    if (index != -1)
                    {
                        // 在参数名数组中，将该参数名对应值写入parameter.txt中
                        fout.WriteLine(x[index]);
                    }
                    else
                    {
                        // 不在参数名数组中，即不是待率定参数，直接写入模型输入文件即可
                        fout.WriteLine(parameter);
                    }
                }
            }
        }

        public void RunModel()
        {
            XAJ(filePath);
        }

        public double PostProcessing()
        {
            simulatedValues = ReadValues(filePath + "Q.txt");

            double NSE = CalculateNSE(simulatedValues, measuredValues);

            return 1 - NSE;
        }

        public List<double> ReadValues(string fileName)
        {
            // 初始化 I/O 变量
            using (var fin = new System.IO.StreamReader(fileName))
            {
                var values = new List<double>();
                string line;
                while ((line = fin.ReadLine()) != null)
                {
                    if (double.TryParse(line, out double val))
                    {
                        values.Add(val);
                    }
                }
                return values;
            }
        }

        public double CalculateNSE(List<double> simulatedValues, List<double> measuredValues)
        {
            double measuredValuesSum = measuredValues.Sum();
            double measuredValuesAvg = measuredValuesSum / measuredValues.Count;

            double numerator = 0.0;
            double denominator = 0.0;

            foreach (double val in measuredValues)
            {
                denominator += Math.Pow(val - measuredValuesAvg, 2);
            }

            for (int index = 0; index < measuredValues.Count; ++index)
            {
                numerator += Math.Pow(simulatedValues[index] - measuredValues[index], 2);
            }

            double NSE = 1 - numerator / denominator;

            return NSE;
        }

        // 构造函数
        public SCEUA()
        {
            m_nopt = 5;
            m_xname = new List<string> { "x1", "x2", "x3", "x4", "x5" };
            m_a = new List<double> { 0, 0, 0, 0, 0 };
            m_bl = new List<double> { -2, -2, -2, -2, -2 };
            m_bu = new List<double> { 2, 2, 2, 2, 2 };

            m_ngs = 5;
            m_npg = 2 * m_nopt + 1;
            m_nps = m_nopt + 1;
            m_alpha = 1;
            m_beta = 2 * m_nopt + 1;
            m_npt = m_ngs * m_npg;

            m_maxn = 3000;
            m_kstop = 10;
            m_pcento = 0.1;
            m_peps = 0.001;

            m_iniflg = true;
            m_iprint = false;
            m_ideflt = 0;

            filePath = "";
        }

        // 析构函数
        ~SCEUA()
        {
        }

        public SCEUA(int argc_, string[] argv_) : this()
        {
            argc = argc_;

            foreach (string arg in argv_)
            {
                argv.Add(arg);
            }

            //string Path = argv[0]; // 获取可执行文件所在路径

            //filePath = Path.Substring(0, Path.LastIndexOf("\\") + 1); // 获取可执行文件所在目录
            filePath = @".//";
        }

        //=============================================================================================

        public void XAJ(string path)
        {
            Parameter parameter = new Parameter(); // 流域参数
            parameter.ReadFromFile(path); // 读取参数文件

            IO io = new IO(); // 流域输入与输出
            io.ReadFromFile(path); // 从文本中导入降雨和蒸发数据，输出流量过程到文本中

            Watershed watershed = new Watershed(); // 流域分单元
            watershed.ReadFromFile(path); // 读取流域分块信息
            watershed.Calculate(io); // 计算各单元流域各时段降雨和水面蒸发量

            State[] state = new State[watershed.GetnW()]; // 每个单元流域的状态
            for (int nw = 0; nw < watershed.GetnW(); nw++)
            {
                state[nw] = new State();
                state[nw].ReadFromFile(path); // 读取模型计算时段长，h
            }

            Evapotranspiration evapotranspiration = new Evapotranspiration(); // 流域蒸散发
            evapotranspiration.SetParameter(parameter); // 设置参数

            Runoff runoff = new Runoff(); // 流域产流
            runoff.SetParameter(parameter); // 设置参数

            Source source = new Source(); // 流域分水源
            source.SetParameter(parameter); // 设置参数

            Confluence confluence = new Confluence(); // 流域汇流
            confluence.SetParameter(parameter); // 设置参数

            Muskingum muskingum = new Muskingum(); // 流域河道汇流
            muskingum.SetParameter(parameter); // 设置参数

            // 逐时段逐单元流域计算
            int nT = io.nrows; // nT表示时段数
            int nW = watershed.GetnW(); // nW为单元流域个数

            io.m_Q = new double[nT]; // 创建数组存储流域出口断面流量过程

            for (int nt = 0; nt < nT; nt++)            
            {
                state[0].m_Q = 0.0; // m_Q置为0

                for (int nw = 0; nw < nW; nw++)
                {
                    state[nw].SetInput(nt, nw, watershed);

                    evapotranspiration.SetState(state[nw]);
                    evapotranspiration.Calculate();
                    evapotranspiration.UpdateState(state[nw]);

                    runoff.SetState(state[nw]);
                    runoff.Calculate();
                    runoff.UpdateState(state[nw]);

                    source.SetState(state[nw]);
                    source.Calculate();
                    source.UpdateState(state[nw]);

                    confluence.SetState(state[nw]);
                    confluence.Calculate();
                    confluence.UpdateState(state[nw]);

                    muskingum.SetState(state[nw]);
                    muskingum.Calculate();
                    muskingum.UpdateState(state[nw]);

                    state[0].m_Q += state[nw].m_O2;
                }

                io.m_Q[nt] = state[0].m_Q;
            }

            // 输出流域出口断面流量过程到文本Q.txt中
            io.WriteToFile(path);
        }
    }
}